Mouse muscle stem cell engraftment: Bulk RNA-seq gene differential analysis

Génome Québec Innovation Centre, McGill University, Montreal, Canada

Prepare analysis workflow

# knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())
# options(
#   readr.show_progress = FALSE,
#   digits = 2,
#   scipen = 8,
#   future.globals.maxSize = +Inf
# )
knitr::opts_knit$set(root.dir = getwd())

library(tidyverse)
library(tximport)
library(DESeq2)
library(vsn)
library(cowplot)
library(glmpca)
library(apeglm)
library(rhdf5)
library(EnhancedVolcano)
library(pheatmap)
library(ensembldb)
library(scater)
library(wesanderson)
library(PoiClaClu)
library(RColorBrewer)
theme_set(theme_bw())

Auxiliary functions

get_deseq_data <- function(target_tissue, meta_dt, data_dir, tx2knownGene) {
  samples_dir <- file.path(
    data_dir,
    meta_dt %>%
      pull(name)
  )
  cdna_dir <- file.path(samples_dir)
  cdna_files <- file.path(cdna_dir, "abundance.h5")
  names(cdna_files) <-
    meta_dt %>%
    pull(sample)


  txi_scaledTPM <- tximport(cdna_files,
    type = "kallisto",
    tx2gene = tx2knownGene,
    countsFromAbundance = c("scaledTPM")
  )

  txi <- tximport(cdna_files,
    type = "kallisto",
    tx2gene = tx2knownGene,
    countsFromAbundance = c("no")
  )
  scaled_tpm <-
    txi_scaledTPM$counts %>%
    as_tibble(rownames = "gene_id")

  gene_counts <-
    txi$counts %>%
    as_tibble(rownames = "gene_id")

  write_tsv(scaled_tpm, paste0(target_tissue, "_engraftment_scaled_tpm.tsv"))
  write_tsv(gene_counts, paste0(target_tissue, "_engraftment_counts.tsv"))

  design_matrix <-
    model.matrix(
      ~ time * age * batch,
      meta_dt %>%
        column_to_rownames("sample")
    )
  colnames(design_matrix) <- c(
    "intercept",
    "engraf",
    "age",
    "intercept_batch",
    "niche",
    "engraf_batch",
    "age_batch",
    "niche_batch"
  )

  design_matrix <- design_matrix[, c(1, 2, 3, 5, 4, 6, 7, 8)]

  dds <-
    DESeqDataSetFromTximport(
      txi,
      meta_dt %>%
        column_to_rownames("sample"),
      design = design_matrix
    )

  return(dds)
}

# We use apeglm shrinkage to preserve the size of large LFC and compute s-values (the estimated rate of false sign). Shrinkage is useful for ranking and visualization, without the need for arbitrary filters on low count genes. https://doi.org/10.1093/bioinformatics/bty895

# The local false sign rate FSR is defined as the posterior probability that the posterior mode (MAP) of beta (log2FC) is of a biologically significant effect size (i.e. abs(beta) > lfcThreshold )

get_lfc <- function(dds, coef_name, lfc_threshold, svalue_threshold) {
  lfc <-
    lfcShrink(dds,
      coef = coef_name,
      type = "apeglm",
      svalue = T,
      lfcThreshold = lfc_threshold
    )

  plotMA(lfc, ylim = c(-22, 8))
  abline(h = c(-1, 1), col = "dodgerblue", lwd = 2)

  lfc <-
    lfc %>%
    as_tibble() %>%
    dplyr::select(-baseMean) %>%
    mutate(
      keep = ((abs(log2FoldChange) > lfc_threshold) & (svalue < svalue_threshold))
    ) %>%
    set_names(c(
      paste0(
        c("lfc_", "lfc_se_", "svalue_", "keep_"),
        coef_name
      )
    ))

  return(lfc)
}

get_lfc_estimates <- function(dds, gene_metadata, coef_names, lfc_threshold, svalue_threshold, batch_threshold) {
  lfcs <-
    map_dfc(coef_names,
      get_lfc,
      dds = dds,
      lfc_threshold = lfc_threshold,
      svalue_threshold=svalue_threshold
    )
  
  coef_names_batch <- c(paste0(c("intercept", coef_names), "_batch"),  paste0("SE_", c("intercept", coef_names), "_batch"))

  lfcs <-
    bind_cols(
      mcols(dds)[, c("baseMean", "baseVar", "dispOutlier" , "maxCooks", "intercept", "SE_intercept", c(coef_names, coef_names_batch))] %>%
        as_tibble(rownames = "gene_id") %>%
        dplyr::select(gene_id, everything()),
      lfcs
    )
  
  lfcs <-
  lfcs %>%
  mutate(
    total_batch= abs(intercept_batch) + abs(niche_batch) + abs(age_batch) +abs(engraf_batch),
    keep_niche  = (keep_niche & abs(niche_batch) < lfc_threshold & total_batch < batch_threshold ),
    keep_age = (keep_age & abs(age_batch) < lfc_threshold & total_batch < batch_threshold  ),
    keep_engraf = (keep_engraf & abs(engraf_batch) < lfc_threshold & total_batch < batch_threshold ),
    EAN = paste0(
      as.integer(keep_engraf),
      as.integer(keep_age),
      as.integer(keep_niche)
    ),
    
    
  )

lfcs <-
  lfcs %>%
  dplyr::select(
    gene_id, EAN, baseMean,
    intercept, paste0("lfc_", coef_names),
    paste0("svalue_", coef_names), everything()
  )

lfcs <-
  lfcs %>%
  left_join(., gene_metadata, by = "gene_id") %>%
  dplyr::select(gene_name, everything()) %>%
  mutate(
    gene_name =
      scater::uniquifyFeatureNames(
        ID = gene_id,
        names = gene_name
      )
  )
  return(lfcs)
}


plot_volcano <- function(data, xvar, yvar, title, lfc_threshold, svalue_thresh, figures_dir) {
  gv <-
    EnhancedVolcano(data,
      lab = data %>% pull(gene_name),
      xlab = bquote(~ italic(Moderated) ~ Log[2] ~ FC),
      ylab = bquote(~ -Log[10] ~ italic(svalue)),
      x = xvar,
      y = yvar,
      col = c("grey30", "forestgreen", "red2", "royalblue"),
      pCutoff = svalue_thresh,
      FCcutoff = lfc_threshold,
      #  legend = c("NS", "Log2FC", "svalue", "svalue & Log2FC"),
      legendLabels = c(
        "NS",
        expression(Log[2] ~ FC),
        "svalue",
        expression(s - value ~ and ~ log[2] ~ FC)
      )
    )

  ggsave(file.path(figures_dir, paste0(title, "_volcano.pdf")), gv)
  return(gv)
}

plot_heatmap <- function(target_genes, title, meta_col, data, show_genes = F, figures_dir) {
  df <-
    left_join(target_genes %>%
      dplyr::select(gene_id, gene_name),
    data,
    by = "gene_id"
    ) %>%
    dplyr::select(-gene_id) %>%
    column_to_rownames("gene_name") %>%
    as.matrix()

  gh <-
    pheatmap(df,
      color = colorRampPalette(rev(RColorBrewer::brewer.pal(n = 11, name = "Spectral")))(100),
      cluster_rows = T,
      show_rownames = show_genes,
      treeheight_row = 0,
      fontsize_col = 4,
      fontsize_row = 6,
      angle_col = "45",
      cluster_cols = FALSE,
      annotation_col = meta_col
    )
  ggsave(file.path(figures_dir, paste0(title, "_heatmap_expression.pdf")), gh)

  return(gh)
}

plot_pca <- function(vsd, title, figures_dir) {
  pcaData <-
    plotPCA(vsd,
      intgroup = c("time", "age", "tissue", "batch", "donor"),
      returnData = TRUE
    )
  percentVar <- round(100 * attr(pcaData, "percentVar"))

  gg_plot <-
    ggplot(
      pcaData %>%
        mutate(condition = paste(age, time, sep = "_")),
      aes(x = PC1, y = PC2, color = time, shape = age, label = name)
    ) +
    geom_point(size = 4) +
    geom_line(aes(group = donor), colour = "grey") +
    #  coord_fixed() +
    geom_text(check_overlap = T, angle = 0, size = 3, vjust = 3, nudge_y = 1.5) +
    xlab(paste0("PC1: ", percentVar[1], "% variance")) +
    ylab(paste0("PC2: ", percentVar[2], "% variance")) +
    coord_fixed() +
    ggtitle("PCA with VST data")
  ggsave(file.path(figures_dir, paste0(title, "_pca_bulk_rnaseq.pdf")), gg_plot)
  gg_plot
}

plot_gene <- function(gene, data, figures_dir) {
  df_plot <-
    data %>%
    dplyr::filter(gene_id == gene)

  gene_name <- df_plot$gene_name[1]
  p <- ggplot(
    df_plot,
    aes(x = time, y = log2_scaled_tpm)
  ) +
    geom_point(aes(color = donor)) +
    geom_line(aes(group = donor)) +
    labs(title = gene_name) +
    facet_grid(batch ~ age)
  ggsave(file.path(figures_dir, paste0(gene_name, "_expression_log2_scaled_tpm.pdf")), p)
  return(p)
}

Set filepaths

data_dir <- "./data/kallisto_out"
gene_metadata_filepath <- "./data/metadata/mm10_ens97_gene_meta.txt"
metadata_filepath <- "./data/metadata/samples_metadata_leg.tsv"

Get tx id versions

library(AnnotationHub)

Loading required package: BiocFileCache
Loading required package: dbplyr

Attaching package: 'dbplyr'

The following objects are masked from 'package:dplyr':

    ident, sql

Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio

Attaching package: 'AnnotationHub'

The following object is masked from 'package:Biobase':

    cache

ah <- AnnotationHub()

DEPRECATION: As of AnnotationHub (>2.23.2), default caching location has changed.
  Problematic cache: /home/rick/.cache/AnnotationHub
  See https://bioconductor.org/packages/devel/bioc/vignettes/AnnotationHub/inst/doc/TroubleshootingTheCache.html#default-caching-location-update
snapshotDate(): 2021-05-18

meta_dt <-
  read_tsv(file.path(metadata_filepath)) %>%
  # arrange(tissue, age, time) %>%
  mutate(
    condition = paste(age, time, sep = "_"),
    age = factor(age, levels = c("young", "aged")),
    time = factor(time, levels = c("T0", "T21")),
    batch = factor(batch, levels = c(1, 2))
  )

Rows: 23 Columns: 8── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (7): name, sample, time, tissue, age, donor, type
dbl (1): batch
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

meta_dt

gene_metadata <-
  read_csv(gene_metadata_filepath) %>%
  dplyr::select(-gene_id_version)

Rows: 55487 Columns: 8── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (6): gene_name, gene_id, seqnames, gene_biotype, gene_id_version, description
dbl (2): gene_len, perc_gene_gc
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

gene_metadata <-
  left_join(tx_meta %>% as_tibble() %>%
    group_by(gene_id) %>%
    dplyr::count(gene_name) %>%
    dplyr::select(-n),
  gene_metadata %>%
    dplyr::select(-gene_name),
  by = "gene_id"
  ) %>%
  mutate(
    gene_name =
      scater::uniquifyFeatureNames(
        ID = gene_id,
        names = gene_name
      )
  )

Load data

dds_leg <- get_deseq_data("leg", meta_dt, data_dir, tx2knownGene)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 
summarizing abundance
summarizing counts
summarizing length
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 
summarizing abundance
summarizing counts
summarizing length
using counts and average transcript lengths from tximport

dds_leg

class: DESeqDataSet 
dim: 36558 23 
metadata(1): version
assays(2): counts avgTxLength
rownames(36558): ENSMUSG00000000001 ENSMUSG00000000003 ... ENSMUSG00000118487 ENSMUSG00000118488
rowData names(0):
colnames(23): yng_t0_01 yng_t0_02 ... aged_t21_05 aged_t21_06
colData names(8): name time ... batch condition

summary_dt_leg <-
  colSums(assay(dds_leg)) %>%
  tibble::enframe(name = "sample", value = "counts")
summary_dt_leg

scaled tpm

scaled_tpm_leg <- read_tsv("leg_engraftment_scaled_tpm.tsv")

Rows: 36558 Columns: 24── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr  (1): gene_id
dbl (23): yng_t0_01, yng_t0_02, yng_t0_03, yng_t0_04, yng_t0_05, yng_t0_06, yng_t21_01, yng_t21_02, yng_t21_04, yng_t21_05, yng_t21_06, aged_t0_01, aged_t0_02, aged_t0_03, aged_t0_04, aged_t0_05, aged_t0_06, aged_t21_01, aged_t21_02, a...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

counts_leg <- read_tsv("leg_engraftment_counts.tsv")

Rows: 36558 Columns: 24── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr  (1): gene_id
dbl (23): yng_t0_01, yng_t0_02, yng_t0_03, yng_t0_04, yng_t0_05, yng_t0_06, yng_t21_01, yng_t21_02, yng_t21_04, yng_t21_05, yng_t21_06, aged_t0_01, aged_t0_02, aged_t0_03, aged_t0_04, aged_t0_05, aged_t0_06, aged_t21_01, aged_t21_02, a...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

scaled_tpm_leg_long <-
  scaled_tpm_leg %>%
  pivot_longer(cols = yng_t0_01:aged_t21_06, names_to = "sample", values_to = "scaled_tpm") %>%
  left_join(., meta_dt %>% dplyr::select(sample, time, age, donor, batch, condition), by = "sample")

scaled_tpm_leg_long <-
  left_join(scaled_tpm_leg_long,
    gene_metadata %>%
      dplyr::select(gene_name, gene_id),
    by = "gene_id"
  ) %>%
  mutate(log2_scaled_tpm = log1p(scaled_tpm))

summary_tpm_leg_long <-
  scaled_tpm_leg_long %>%
  group_by(gene_id, condition) %>%
  summarize(
    mean_tpm = mean(log2_scaled_tpm),
    sd_tpm = sd(log2_scaled_tpm),
    frac_zeros = mean(scaled_tpm == 0)
  ) %>%
  mutate(cv = sd_tpm / mean_tpm) %>%
  group_by(gene_id) %>%
  mutate(
    low_frac_zeros = sum(frac_zeros <= 0.2),
    sd_group_tpm = sd(mean_tpm),
    mean_group_tpm = mean(mean_tpm),
    min_tpm = min(mean_tpm),
    max_tpm = max(mean_tpm)
  ) %>%
  dplyr::filter(low_frac_zeros > 0) %>%
  # dplyr::filter(min_tpm >0.1) %>%
  dplyr::filter(max_tpm > log(30))

`summarise()` has grouped output by 'gene_id'. You can override using the `.groups` argument.

summary_tpm_leg_long

ggplot(summary_tpm_leg_long %>%
  group_by(gene_id) %>%
  dplyr::slice(1)) +
  geom_point(aes(x = min_tpm, y = max_tpm - min_tpm), size = 0.5)

target_genes <-
  summary_tpm_leg_long %>%
  pull(gene_id) %>%
  unique()

keep <- rownames(counts(dds_leg)) %in% target_genes
dds_leg <- dds_leg[keep, ]

Generalized Principal Component Analysis

set.seed(11678)
gpca <- glmpca(scaled_tpm_leg %>%
  dplyr::filter(gene_id %in% target_genes) %>%
  column_to_rownames("gene_id"),
L = 2,
optimizer = c("fisher"),
fam = c("poi"),
X = model.matrix(~ -1 + batch, meta_dt %>% dplyr::select(batch))
)
gpca.dat <- gpca$factors
gpca.dat <- cbind(gpca.dat, meta_dt[, c("age", "time", "batch", "donor")])
gpca.dat

gg_pca <-
  ggplot(
    gpca.dat %>% as_tibble(rownames = "sample") %>%
      mutate(condition = paste(age, time, sep = "_")),
    aes(x = dim1, y = dim2, color = condition, shape = batch, label = sample)
  ) +
  geom_point(size = 4) +
  geom_line(aes(group = donor), colour = "grey") +
  #  coord_fixed() +
  geom_text(check_overlap = T, angle = 0, size = 3, vjust = 2, nudge_y = 1.5) +
  # facet_grid(~tissue) +
  ggtitle("glmpca - Generalized PCA")
gg_pca

colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
poisd <- PoissonDistance(t(counts(dds_leg)))
samplePoisDistMatrix <- as.matrix( poisd$dd )
rownames(samplePoisDistMatrix) <- colnames(counts(dds_leg))
colnames(samplePoisDistMatrix) <- colnames(counts(dds_leg))
pheatmap(samplePoisDistMatrix,
         clustering_distance_rows = poisd$dd,
         clustering_distance_cols = poisd$dd,
         col = colors)

vsd <- vst(dds_leg, blind = T)

using 'avgTxLength' from assays(dds), correcting for library size

pcaData <- plotPCA(vsd, intgroup = c("condition", "batch"), returnData = TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))

ggplot(pcaData, aes(x = PC1, y = PC2, color = condition, shape = batch)) +
  geom_point(size =3) +
  xlab(paste0("PC1: ", percentVar[1], "% variance")) +
  ylab(paste0("PC2: ", percentVar[2], "% variance")) +
  coord_fixed() +
  ggtitle("PCA with VST data")

design_matrix <-
  model.matrix(
    ~1,
    meta_dt %>%
      column_to_rownames("sample")
  )
colnames(design_matrix) <- c("intercept")
dds_leg <- DESeq(dds_leg, test = "LRT", reduced = design_matrix)

using supplied model matrix
estimating size factors
using 'avgTxLength' from assays(dds), correcting for library size
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

plotDispEsts(dds_leg)

res <- results(dds_leg)
W <- res$stat
maxCooks <- apply(assays(dds_leg)[["cooks"]],1,max)
idx <- !is.na(W)
plot(rank(W[idx]), maxCooks[idx], xlab="rank of Wald statistic", 
     ylab="maximum Cook's distance per gene",
     ylim=c(0,5), cex=.4, col=rgb(0,0,0,.3))
m <- ncol(dds_leg)
p <- 3
abline(h=qf(.99, p, m - p))

plot(res$baseMean+1, -log10(res$pvalue),
     log="x", xlab="mean of normalized counts",
     ylab=expression(-log[10](pvalue)),
     ylim=c(0,30),
     cex=.4, col=rgb(0,0,0,.3))

lfc_threshold <- 1
svalue_thresh <- 0.05
batch_threshold <- 4
coef_names <- c("engraf", "age", "niche")
lfcs <- 
  get_lfc_estimates(dds_leg,
                    gene_metadata,
                    coef_names,
                    lfc_threshold,
                    svalue_thresh,
                    batch_threshold)

using 'apeglm' for LFC shrinkage. If used in published research, please cite:
    Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
    sequence count data: removing the noise and preserving large differences.
    Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
computing FSOS 'false sign or small' s-values (T=1)
thresholding s-values on alpha=0.005 to color points

using 'apeglm' for LFC shrinkage. If used in published research, please cite:
    Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
    sequence count data: removing the noise and preserving large differences.
    Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
computing FSOS 'false sign or small' s-values (T=1)
thresholding s-values on alpha=0.005 to color points

lfcs

Plot

p1 <- ggplot(
  lfcs,
  aes(x = intercept_batch, y = intercept)
) +
  geom_point(size = .7, alpha = 0.3) +
  geom_vline(xintercept = c(-lfc_threshold, lfc_threshold)) +
  geom_hline(yintercept = c(-1, 1))
# ggsave(file.path(figures_dir, "niche_vs_age_lfc_leg.pdf"), p1)
p1

p2 <- ggplot(
  lfcs,
  aes(x = engraf_batch, y = lfc_engraf, color = svalue_engraf < 0.1)
) +
  geom_point(size = .7, alpha = 0.3) +
  geom_vline(xintercept = c(-lfc_threshold , lfc_threshold)) +
  geom_hline(yintercept = c(-1, 1))
# ggsave(file.path(figures_dir, "niche_vs_age_lfc_leg.pdf"), p1)
p2

p3 <- ggplot(
  lfcs,
  aes(x = age_batch, y = lfc_age, color = svalue_age < 0.10)
) +
  geom_point(size = .7, alpha = 0.3) +
  geom_vline(xintercept = c(-lfc_threshold , lfc_threshold )) +
  geom_hline(yintercept = c(-1, 1))
# ggsave(file.path(figures_dir, "niche_vs_age_lfc_leg.pdf"), p1)
p3

p4 <- ggplot(
  lfcs,
  aes(x = niche_batch, y = lfc_niche, color = svalue_niche < 0.05)
) +
  geom_point(size = .7, alpha = 0.3) +
  geom_vline(xintercept = c(-lfc_threshold, lfc_threshold )) +
  geom_hline(yintercept = c(-1, 1))
# ggsave(file.path(figures_dir, "niche_vs_age_lfc_leg.pdf"), p1)
p4

table(lfcs$EAN)

  000   001   010   011   100   101   110   111 
15204    30    59    12    46     2     6     3 

 ggplot(
  lfcs,
  aes(x = total_batch)
) +
  geom_histogram(bins=60) +
  scale_x_log10() +
  geom_vline(xintercept = c(3))

Write to file

write_tsv(lfcs, "leg_two_batches_moderated_lfc_1_svalue_05_marked.txt")

lfcs_hits <-
  lfcs %>%
  dplyr::filter(lfcs$EAN != "000")
lfcs_hits

Hits by group

target_genes_leg_001 <-
  lfcs_hits %>%
  dplyr::filter(EAN %in% c("001")) %>%
  arrange(svalue_niche)
target_genes_leg_001

figures_dir <- "./figures/hits/leg_001"
pp1 <- map(target_genes_leg_001 %>% pull(gene_id), plot_gene, data = scaled_tpm_leg_long, figures_dir)

Saving 7 x 7 in image

target_genes_leg_011 <-
  lfcs_hits %>%
  dplyr::filter(EAN %in% c("011")) %>%
  arrange(-baseMean)
target_genes_leg_011

figures_dir <- "./figures/hits/leg_011"
pp2 <- map(target_genes_leg_011 %>% pull(gene_id), plot_gene, data = scaled_tpm_leg_long, figures_dir)

Saving 7 x 7 in image

target_genes_leg_101 <-
  lfcs_hits %>%
  dplyr::filter(EAN %in% c("101"))
target_genes_leg_101

figures_dir <- "./figures/hits/leg_101"
pp3 <- map(target_genes_leg_101 %>% pull(gene_id), plot_gene, data = scaled_tpm_leg_long, figures_dir)

Saving 7 x 7 in image

target_genes_leg_111 <-
  lfcs_hits %>%
  dplyr::filter(EAN %in% c("111"))
target_genes_leg_111

figures_dir <- "./figures/hits/leg_111"
pp4 <- map(target_genes_leg_111 %>% pull(gene_id), plot_gene, data = scaled_tpm_leg_long, figures_dir)

Saving 7 x 7 in image

target_genes_leg_100 <-
  lfcs_hits %>%
  dplyr::filter(EAN %in% c("100"))
target_genes_leg_100

figures_dir <- "./figures/hits/leg_100"
pp5 <- map(target_genes_leg_100 %>% pull(gene_id), plot_gene, data = scaled_tpm_leg_long, figures_dir)

Saving 7 x 7 in image

target_genes_leg_010 <-
  lfcs_hits %>%
  dplyr::filter(EAN %in% c("010"))
target_genes_leg_010

figures_dir <- "./figures/hits/leg_010"
pp6 <- map(target_genes_leg_010 %>% pull(gene_id), plot_gene, data = scaled_tpm_leg_long, figures_dir)

Saving 7 x 7 in image

target_genes_leg_110 <-
  lfcs_hits %>%
  dplyr::filter(EAN %in% c("110"))
target_genes_leg_110

figures_dir <- "./figures/hits/leg_110"
pp7 <- map(target_genes_leg_110 %>% pull(gene_id), plot_gene, data = scaled_tpm_leg_long, figures_dir)

Saving 7 x 7 in image

genes_to_plot <- c("Pax7", "Myf5", "Myod1")
figures_dir <- "./figures/additional"
map(lfcs %>%
  dplyr::filter(gene_name %in% genes_to_plot) %>% pull(gene_id), 
  plot_gene,
  data = scaled_tpm_leg_long,
  figures_dir)

Saving 7 x 7 in image

[[1]]

[[2]]

[[3]]

meta_df_leg <-
  as.data.frame(colData(dds_leg)[, c("age", "time", "batch")])

dt_plot_leg <-
  scaled_tpm_leg %>%
  mutate_if(is.numeric, ~ log1p(.))

figures_dir <- "./figures/hits"
plot_heatmap(lfcs_hits, title = "log2_scaled_tpm_leg_hits", meta_df_leg, dt_plot_leg, show_genes = F, figures_dir)

figures_dir <- "./figures/hits"
plot_heatmap(target_genes_leg_001, title = "log2_scaled_tpm_leg_hits_001", meta_df_leg, dt_plot_leg, show_genes = T, figures_dir)

figures_dir <- "./figures/hits"
plot_heatmap(target_genes_leg_100, title = "log2_scaled_tpm_leg_hits_100", meta_df_leg, dt_plot_leg, show_genes = T, figures_dir)

figures_dir <- "./figures/hits"
plot_heatmap(target_genes_leg_010, title = "log2_scaled_tpm_leg_hits_010", meta_df_leg, dt_plot_leg, show_genes = T, figures_dir)

svalue_thresh <- 0.05
figures_dir <- "./figures/volcanos"
plot_volcano(lfcs, "lfc_niche", "svalue_niche", "niche", lfc_threshold, svalue_thresh, figures_dir)

One or more p-values is 0. Converting to 10^-1 * current lowest non-zero p-value...Ignoring unknown parameters: xlim, ylimSaving 7 x 7 in image

plot_volcano(lfcs, "lfc_age", "svalue_age", "age", lfc_threshold, svalue_thresh , figures_dir)

Ignoring unknown parameters: xlim, ylimSaving 7 x 7 in image

plot_volcano(lfcs, "lfc_engraf", "svalue_engraf", "engraf", lfc_threshold, svalue_thresh, figures_dir)

Ignoring unknown parameters: xlim, ylimSaving 7 x 7 in image

sessionInfo()

R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] AnnotationHub_3.0.2         BiocFileCache_2.0.0         dbplyr_2.1.1                RColorBrewer_1.1-2          PoiClaClu_1.0.2.1           wesanderson_0.3.6           scater_1.20.1               scuttle_1.2.1              
 [9] SingleCellExperiment_1.14.1 ensembldb_2.16.4            AnnotationFilter_1.16.0     GenomicFeatures_1.44.2      AnnotationDbi_1.54.1        pheatmap_1.0.12             EnhancedVolcano_1.10.0      ggrepel_0.9.1              
[17] rhdf5_2.36.0                apeglm_1.14.0               glmpca_0.2.0                cowplot_1.1.1               vsn_3.60.0                  DESeq2_1.32.0               SummarizedExperiment_1.22.0 Biobase_2.52.0             
[25] MatrixGenerics_1.4.3        matrixStats_0.61.0          GenomicRanges_1.44.0        GenomeInfoDb_1.28.4         IRanges_2.26.0              S4Vectors_0.30.2            BiocGenerics_0.38.0         tximport_1.20.0            
[33] forcats_0.5.1               stringr_1.4.0               dplyr_1.0.7                 purrr_0.3.4                 readr_2.1.2                 tidyr_1.1.4                 tibble_3.1.6                ggplot2_3.3.5              
[41] tidyverse_1.3.1            

loaded via a namespace (and not attached):
  [1] utf8_1.2.2                    tidyselect_1.1.1              RSQLite_2.2.9                 grid_4.1.2                    BiocParallel_1.26.2           munsell_0.5.0                 ScaledMatrix_1.0.0           
  [8] ragg_1.2.1                    preprocessCore_1.54.0         withr_2.4.3                   colorspace_2.0-2              filelock_1.0.2                ggalt_0.4.0                   knitr_1.37                   
 [15] rstudioapi_0.13               Rttf2pt1_1.3.9                labeling_0.4.2                bbmle_1.0.24                  GenomeInfoDbData_1.2.6        farver_2.1.0                  bit64_4.0.5                  
 [22] coda_0.19-4                   vctrs_0.3.8                   generics_0.1.1                xfun_0.29                     R6_2.5.1                      ggbeeswarm_0.6.0              rsvd_1.0.5                   
 [29] locfit_1.5-9.4                bitops_1.0-7                  rhdf5filters_1.4.0            cachem_1.0.6                  DelayedArray_0.18.0           assertthat_0.2.1              vroom_1.5.7                  
 [36] promises_1.2.0.1              BiocIO_1.2.0                  scales_1.1.1                  beeswarm_0.4.0                gtable_0.3.0                  beachmat_2.8.1                ash_1.0-15                   
 [43] affy_1.70.0                   rlang_1.0.0                   genefilter_1.74.1             systemfonts_1.0.3             splines_4.1.2                 rtracklayer_1.52.1            extrafontdb_1.0              
 [50] lazyeval_0.2.2                broom_0.7.12                  BiocManager_1.30.16           yaml_2.2.2                    modelr_0.1.8                  backports_1.4.1               httpuv_1.6.5                 
 [57] extrafont_0.17                tools_4.1.2                   affyio_1.62.0                 ellipsis_0.3.2                jquerylib_0.1.4               Rcpp_1.0.8                    plyr_1.8.6                   
 [64] sparseMatrixStats_1.4.2       progress_1.2.2                zlibbioc_1.38.0               RCurl_1.98-1.5                prettyunits_1.1.1             viridis_0.6.2                 haven_2.4.3                  
 [71] fs_1.5.2                      magrittr_2.0.2                reprex_2.0.1                  mvtnorm_1.1-3                 ProtGenerics_1.24.0           evaluate_0.14                 mime_0.12                    
 [78] hms_1.1.1                     xtable_1.8-4                  XML_3.99-0.8                  emdbook_1.3.12                readxl_1.3.1                  gridExtra_2.3                 compiler_4.1.2               
 [85] biomaRt_2.48.3                bdsmatrix_1.3-4               maps_3.4.0                    KernSmooth_2.23-20            crayon_1.4.2                  htmltools_0.5.2               later_1.3.0                  
 [92] tzdb_0.2.0                    geneplotter_1.70.0            lubridate_1.8.0               DBI_1.1.2                     proj4_1.0-11                  MASS_7.3-55                   rappdirs_0.3.3               
 [99] Matrix_1.4-0                  cli_3.1.1                     pkgconfig_2.0.3               GenomicAlignments_1.28.0      numDeriv_2016.8-1.1           xml2_1.3.3                    annotate_1.70.0              
[106] bslib_0.3.1                   vipor_0.4.5                   XVector_0.32.0                rvest_1.0.2                   digest_0.6.29                 Biostrings_2.60.2             rmarkdown_2.11               
[113] cellranger_1.1.0              DelayedMatrixStats_1.14.3     restfulr_0.0.13               curl_4.3.2                    shiny_1.7.1                   Rsamtools_2.8.0               rjson_0.2.21                 
[120] lifecycle_1.0.1               jsonlite_1.7.3                Rhdf5lib_1.14.2               BiocNeighbors_1.10.0          viridisLite_0.4.0             limma_3.48.3                  fansi_1.0.2                  
[127] pillar_1.6.5                  lattice_0.20-45               ggrastr_1.0.1                 KEGGREST_1.32.0               fastmap_1.1.0                 httr_1.4.2                    survival_3.2-13              
[134] interactiveDisplayBase_1.30.0 glue_1.6.1                    png_0.1-7                     BiocVersion_3.13.1            bit_4.0.4                     sass_0.4.0                    stringi_1.7.6                
[141] blob_1.2.2                    textshaping_0.3.6             BiocSingular_1.8.1            memoise_2.0.1                 irlba_2.3.5